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A  CLASS  OF  MONOTONIC  SHOCK-CAPTURING  DIFFERENCE  SCHEMES 
0.  Introduction 

A  broad  class  of  gas-dynamic  problems  can  be  treated  within  the 
framework  of  ideal  gas  theory,  i.e.,  assuming  the  gas  is  inviscid  and 
thermally  nonconductive.  The  differential  equations  of  gas  dynamics 
follow  from  the  laws  of  conservation  of  mass,  momentum,  and  energy 
in  integral  form  assuming  continuous  differentiability  of  the  fluid 
variables,  and  constitute  a  system  of  quasi-linear  hyperbolic  equations. 

It  is  well  known  that  in  general  these  may  posess  discontinuous  solutions 
even  when  smooth  initial  data  are  prescribed*.  Physically,  the  presence 
of  discontinuities  in  the  solution  typically  signals  the  appearance  of 
a  shock  wave.  There  are  two  possibilities  for  correctly  describing 
discontinuous  flows  within  the  framework  of  ideal  gas  theory.  The 
first  consists  of  breaking  up  all  the  regions  In  the  problem  into 
subregions  of  smooth  flow,  which  are  described  by  the  differential 
equations  of  gas  dynamics,  while  the  discontinuities  (boundaries  of  the 
subregions),  are  described  by  conservation  conditions.  It  is  important- 

to  distinguish  weak  discontinuities  (discontinuities  in  derivatives  of 
the  fluid  variables),  tangential  discontinuties,  and  finally,  strong 

discontinuities  (shock  waves) .  The  second  approach  consists  in  utilizing 

the  conservation  conditions  in  integral  form,  which  allows  discontinuous 
solutions.  For  historical  reasons  (the  differential  equations  were 

derived  earlier),  these  are  referred  to  as  "generalized,*1  iu  contrast 
to  the  continuous  classical  equations.  Ve  note  that  consideration  of 

generalized  solutions  extends  the  class  of  possible  solutions,  and  it 
is  therefore  necessary  to  make  use  of  additional  considerations  in  order 

to  determine  the  correct  solution.  For  example  the  requirement  that 
entropy  not  decrease  at  a  discontinuity  permits  us  to  exclude 
Manuscript  submitted  May  12,  1981. 
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rarefaction  shocks  in  the  flow  of  a  perfect  gas  (from  a  formal  point 

of  view  such  a  discontinuity  is  unstable) . 

In  numerical  modelling  of  gas  motions  in  the  presence  of  shock  waves 

both  of  the  above  approaches  are  possible.  In  the  former,  the  physical 

jump  conditions  on  the  discontinuities  dividing  the  regions  of  continuous 

flow  are  used  both  in  order  to  obtain  the  conditions  connecting  the  fluid 

variables  on  the  two  sides  of  the  discontinuity  and  to  determine  the 

motion  of  the  computational  mesh  points  following  the  discontinuity. 

Changes  in  the  topology  of  the  lines  of  discontinuity  can  be  followed 

either  using  a  moving  curvilinear  mesh  or  by  means  of  an  algorithm  based 

on  a  fixed  mesh.  The  first  method  is  most  completely  worked  out  in  the 

2 

paper  by  S.  K.  Godunov  et  al.  while  the  second  has  been  developed  by 
,  3,4 

Moretti  et  al. 

The  shortcoming  of  both  approaches  is  the  inhomogeneity  of  the 
resulting  difference  schemes  and  consequent  complicated  structure  of 
numerical  algorithm.  Both  approaches,  as  a  rule,  demand  an  a  priori  know¬ 
ledge  of  the  flow  pattern  during  the  calculation. 

When  it  is  impossible  to  know  in  advance  the  flow  pattern  ,dr  when 
it  is  changing  qualitatively  with  time,  it  is  more  convenient  to  make  use 
of  "shock  capturing"  difference  schemes,  amounting  to  difference  approxi¬ 
mations  to  the  integral  form  of  the  conservation  laws  for  every  computa¬ 
tional  cell.  The  form  of  the  difference  equations  does  not  depend  on 
the  character  of  the  flow  or  the  position  of  the  possible  shocks,  and 
therefore  such  schemes  are  homogeneous. [Apparently  "homogeneous" 

•c 

means  universal,  i.e.,  not  problem-dependent. — DLB]  In  this  approach 
shocks  appear  as  regions  of  abrupt  variation  of  the  fluid  quantities. 
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Shock-capturing  schemes  can  also  be  based  on  the  differential  equations. 
For  this  purpose  small  corrections  are  introduced  in  the  equations, 
typically  of  nonlinear  form,  analogous  to  physical  viscosity.  The 
equations  assume  parabolic  form  and  permit  a  smooth  solution  which 
tends  to  the  solution  of  the  original  system  as  the  "artificial"  vis¬ 
cosity  vanishes.  However,  difference  analogs  of  the  conservation  laws 
are  satisfied  in  this  technique,  generally  speaking,  only  approximately. 
In  difference  schemes  based  on  the  integral  form,  the  mass,  momentum, 
and  energy  conservation  laws  are  sati^'-'ed  for  every  computational  cell 
to  roundoff.  Such  schemes  are  describable  in  terms  of  fluxes  across 
boundaries  of  the  computational  cells,  and  so  the  conservation  laws 
for  each  computational  region  are  algebraic  consequences  of  the  con¬ 
servation  laws  for  cells,  i.e.,  the  schemes  are  conservative.^ 

We  note  that  the  first  approach,  as  a  rule,  distinguishes  only 
simple  shocks,  while  the  other  types  of  discontinuity  are  calculated 
as  in  the  shock-capturing  method.  A  difference  scheme  employed  for 
solving  practical  problems  is  required  to  be  accurate,  that  is,  it 
must  be  able  to  describe  flows  on  comparatively  coarse  meashes;  it  must 
be  economical;  and  it  should  be  simple  to  employ.  The  accuracy  of  the 
scheme  for  calculating  smooth  flows  is  determined  by  the  order  of  the 
approximation.  In  the  neighborhood  of  a  shock,  the  change  in  the  fluid 
variables  is  comparable  with  their  magnitude,  so  that  the  concept  of 
the  order  of  the  approximation  becomes  meaningless.  It  is  well  known 
that  schemes  of  the  first  type  lead  to  smearing  of  the  shock  because 
of  the  strong  numerical  viscosity.  Schemes  of  second  and  higher  order 
give  rise  to  significantly  less  smearing  of  shock  runs,  but  are  typically 
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nonmonotonic.  One  of  the  principle  reasons  for  the  appearance  of  un- 
physical  oscillations  in  the  neighborhood  of  a  shock  is  the  nonvanishing 
dispersion  of  the  difference  scheme,  which  leads  to  misrepresentations 
of  the  form  and  speed  of  propagation  of  physical  disturbances,  especially 
at  short  wavelengths.  Nonmonotonicity  may  also  be  produced  by  the  non¬ 
linearity  inherent  in  real  problems.® 

Artificial  viscosity  is  introduced  to  suppress  unphysical  oscil¬ 
lations.  Most  methods  for  achieving  this,  however,  do  not  yield  a 
monotonic  scheme,  and  consequently  make  the  localization  of  shocks  worse 
and  strongly  smear  large  gradients.  This  problem  becomes  especially 
severe  in  calculations  involving  strongly  shocked  flow  arising  from 
nonsteady  interaction  of  shock  waves  with  one  another  and  with  obstacles. ^ 

Recently  a  number  of  methods  using  nonlinear  limiters  of  various  sorts 

7-14 

to  filter  unphysical  oscillations  have  appeared.  Many  of  these  are 

highly  effective,  so  that  the  assertion  that  "good  schemes  of  first  and 

2 

of  higher  order  smear  shock  fronts  in  practically  the  same  degree" 
cannot  be  accepted  as  correct. 

It  is  convenient  to  choose  an  actual  difference  scheme  in  several 
stages.  First  the  dispersion  and  dissipation  properties  are  studied 
in  the  linear  approximation  using  the  method  of  differential  approxi¬ 
mations  or  harmonic  analysis  for  the  simplest  transport  equations.  It 
is  essential  to  test  the  schemes  selected  through  linear  analysis  on 
simple  one-dimensional  problems  having  an  exact  solution.  The  final 
stage  of  development  of  any  difference  scheme  must  be  its  verification 
on  model  two-dimensional  problems. 
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In  the  present  paper  a  nonlinear  conservative  smoothing  procedure 
is  proposed  in  order  to  achieve  monotonicity  in  explicit  schemes  of 
higher  order.'  Although  this  smoothing  procedure  is  applicable  to 
arbitrary  schemes,  the  underlying  transport  algorithm  is  taken  here  to 
be  second-order  HacCormack  differencing.^  The  results  of  the  linear 
analysis  of  the  dispersion  and  dissipation  properties  of  this  and  a 
number  of  other  widely  used  schemes  are  presented  in  Ref.  16.  Below, 
the  second  stage  in  the  choice  of  an  optimal  scheme  is  considered  in 
detail.  This  is  conveniently  divided  into  two  parts.  The  first  is  an 
investigation  of  the  solution  of  the  linear  advectlon  equation,  which 
facilitates  an  understanding  of  the  nonlinear  properties  of  the  scheme. 

The  second  part  is  a  calculation  of  the  flow  resulting  from  the  evolution 
of  a  one-dimensional  gas  dynamic  shock.  In  conclusion,  as  an  illus¬ 
tration  of  the  possibilities  of  the  proposed  numerical  method,  a 
three-dimensional  problem  involving  the  interaction  of  a  plane  shock 
wave  with  an  obstacle  is  discussed. 
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1.  Construction  of  second  order  monotonic  difference  schemes 

An  effective  device  for  constructing  monotonic  difference  schemes 

having  second-order  accuracy  for  smooth  flows  was  presented  by  Boris 
11  12 

and  Book.  '  It  consists  of  two  stages.  In  the  first  a  large  numeri¬ 
cal  diffusion  is  introduced  into  the  solution,  which  guarantees  the 
monotonicity  of  the  scheme.  In  the  second  stage  this  diffusion  is 
canceled  wherever  doing  so  does  not  introduce  new  unphysical  extrema  or 

accentuate  existing  ones.  By  means  of  this  approach  a  family 
of  FCT  (Flux-Corrected  Transport)  methods  was  constructed. 

We  consider  the  simple  example  of  the  Cauchy  problem  for  the 
linear  advection  equation 


where 


3f 

3t 


+  v 


3f 

3x 


0 


(1.1) 


f(0,x)  «  fQ(x). 

The  exact  solution  has  the  form  f(t,x)  =  f^Cx  “  vt) .  Let  f^  be  the 
value  of  the  function  at  the  i^  grid  point  of  a  uniform  mesh  on  the  nth 
time  level;  let  I  +  T  be  the  operator  which  propagates  the  solution 
from  the  nth  to  the  (n+l)st  level.  The  form  of  T  depends  on  the 
particular  difference  scheme  used.  We  define  a  diffusion  operator  D  by 


Dfi  ‘  W  ‘  ’  Q6W  -  Q!Ei-l/2 

where 

6fi+l/2  =  fi+l  “  fr 


(1.2) 
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Then  the  first  (diffusion)  stage  can  be  represented  in  the  fora 

*=  (1+T+D)f°. 

Now  we  cancel  the  diffusion  in  such  a  way  that  in  the  resulting  solu- 
tion  no  new  extrema  in  comparison  with  f ^  appear,  and  those  already 
present  are  not  amplified.  For  this  purpose  in  the  antidiffusion 
operator  A,  given  by 

M i  =  "  (<?i+l/2  “  q?i-l/2>  ’ 

we  limit  the  fluxes  “  ®6f±+l/2  accorditlg  to  the  formula 

*1+1/2  ‘  S  '  0311  {0>  I'Pi+l/zl' 

S  5T1+3/2)>.  (1-3) 

where  S  =  sign  cp^+1/2* 

The  transition  from  the  nth  to  the  (n+l)st  time  level  has  the  fora 

f”+1  =  (1+A)^  =  (1+A)(l+T+D)f^  .  (1.4) 

This  technique — the  explicit  cancellation  of  the  diffusion — leads  to 
retention  of  some  of  the  diffusion  in  smooth  regions  [i.e.,  where  the 
limiter  (1.3)  does  not  operate  and  A  =?  -  D] ,  even  when  T  =  0.  Conse- 

4 

quently  Boris  and  Book  proposed  two  other  algorithms: 
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"phoenical*1 

fj+1  *  [  (1+A)  (1+T)  +D]fJ  (1.5) 

and 

Implicit 

fj+1  =  (l+A)~1(l+T-HD)f°>  (1.6) 

which  permit  complete  cancellation  of  the  diffusion  when  T  =  0  in 

e 

smooth  regions.  We  note  that  using  the  algorithm  (1.4) -(1.6)  leads  to 
substantial  improvement  in  the  results  in  comparison  with  widely  used 
schemes . 

In  the  SHASTA^  scheme  the  coefficient  Q  **  1/8.  Subsequently 

Boris  and  Book  investigated  in  detail  the  influence  of  the  value  of 

magnitude  of  Q  on  the  dissipative  properties  of  this  scheme  and  deter- 

12 

mined  the  optimum  dependence  of  Q  on  the  Courant  number.  However, 
practically  speaking,  this  optimization  is  meaningful  only  for  a  linear 
equation  with  constant  coefficients.  For  the  equations  of  gas-dynamics, 
as  calculations  reveal,  good  results  are  obtained  for  Q  between 
1/10  and  1/6. 

One  can  point  out  three  factors  responsible  for  the  success  of  the 
method  of  Boris  and  Book.  First,  the  diffusion  stage,  as  shown  by 
linear  analysis,  substantially  improves  the  dispersion  properties  of 
the  scheme.  Second,  the  diffusion  is  introduced  in  a  conservative 
manner,  and,  finally,  it  is  done  so  as  Co  permit  excellent  localization 
in  the  neighborhood  of  the  shocks.  Evidently  it  is  precisely  these 
last  two  properties  which  are  most  important.  In  Ref.  13  it  w£s  . 
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proposed  to  interpret  the  method  of  Boris  and  Book  as  a  nonlinear  con¬ 
servative  smoothing  procedure  and  to  introduce  diffusion  at  the  (n+l)st 
time  level  according  to 


f"  -  (1+A)(l+D)(l+T)fj 
or 

f"+1  *=  (1+A+D)(l+T)f° 


(1.7) 


(1.8) 


This  procedure  fails  to  improve  the  dispersive  properties  and  falls 
somewhat  short  in  quality  of  the  results  of  the  algorithms  (1.4)-(1.6). 
However,  in  contrast  to  the  original  method  of  Boris  and  Book,  which  is 
essentially  one-dimensional  in  character  and  which  in  two- 
dimensional  calculations  can  only  be  applied  via  the  method  of  coordi¬ 
nate  timestep-splitting  of  the  the  difference  equations,  smoothing  in 
the  form  of  Eqs.  (1.7)-(1.8)  is  easily  incorporated  in  essentially 
two-dimensional  schemes.  For  this  purpose,  after  each  timestep  successive 
smoothing  operations  are  carried  out  with  respect  to  each  of  the  coordi¬ 
nates.  Such  smoothing  has  turned  out  to  be  effective  in  solving  a 
whole  series  of  problems. ^  ^ 

In  the  present  paper  a  simpler  smoothing  procedure  is  proposed, 
having  a  local  character  while  retaining  the  conservative  property  of 
the  difference  scheme.  A  constant  diffusion  is  introduced  only  in 
regions  which  are  not  monotonic,  i.e.,  where  the  solution  has  a  numerical 
ripple.  As  a  result  of  this,  just  as  in  the  method  of  Boris  and  Book, 
there  is  a  certain  amount  of  smearing  of  physical  extrema,  but  for  the 

4 

most  part  only  oscillations  produced  by  the  numerical  scheme  are 
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removed.  The  conservative  property  is  assured  by  introducing  the 
diffusion  in  the  form  of  fluxes  across  the  cell  boundaries.  Using 
the  notation  introduced  above,  we  can  write  the  smoothing  at  the  nth 
time  level  as 

f"+1  -  (1+T)f*  +  D*fJ,  (1.9) 

and  on  the  (n+l)st  time  level  in  the  form 

l 

f^1  -  (1+D*)(l+T)f^.  (1.10) 

In  these  expressions 

D  fi  "  9i+i/2  “  Cpi-l/2’ 

where 

*  9i+l/2’  6fi+l/2  6fi+3/2  <  0  or  <Sfi+l/2  6fi-l/2  <  0 

Vl/2  = 

0  otherwise. 

In  order  to  code  this  prescription  it  is  convenient  to  search  the  entire- 
array  f recording  the  boundaries  across  which  the  flux  is  nonvanishing , 
and  then  in  a  second  pass  calculate  the  nonzero  fluxes. 

This  procedure  for  introducing  artificial  viscosity  runs  substan~ 
tially  faster  than  the  method  of  Boris  and  Book,  while  yielding  com¬ 
parable  results.  Note  that,  in  contrast  with  the  algorithms  (1.4)-(1.7), 
smoothing  according  to  the  prescription  (lc8)-(1.10)  in  regions  where 
the  fluid  quantities  vary  monotonically  in  general  leaves  the  solution 
unchanged. 


2.  Solution  of  the  linear  transport  equation 

In  this  and  the  following  section  the  proposed  numerical  method 

is  compared  with  a  number  of  different  schemes,  using  a  one  dimensional 

22 

model  problem  having  an  exact  solution.  The  schemes  of  S.  T.  Godunov  , 

15  23 

MacCormack  ,  V.  V.  Eremin  and  Yu.  M.  Lipnitskii  ,  of  first,  second 

and  third  order,  respectively,  are  considered  and  also  the  non-»linear 

12 

scheme  of  Boris  and  Book,  SHASTX  ,  embodying  the  algorithm  (1.5) .  These 
schemes  are  widely  employed  and  have  performed  well  in  practice;  therefore 
the  present  analysis  permits  one  to  carry  out  an  objective  evaluation  of 
the  proposed  method. 

Let  us  approximate  equation  (1.1)  on  a  uniform  mesh  with  Ax  ■  1. 

The  computational  mesh  contains  a  total  of  100  cells.  On  the  boundaries 
periodicity  conditions  are  Imposed.  The  initial  conditions  are  given  in 
the  form 

!0.5,  x  >  L 

<p(x)  ,  1  &  x  &  L 

where  L  =  20  and 

cp(x)  =2  (2.1) 

or 

<p(x)  -  0.5  +  0.075x.  (2.2) 

In  Fig..  1  are  shown  the  results  of  the  calculations  after  800 
time  steps  with  Courant  number  v  =  0.2,  i.e.  for  a  total  running 
time  t  =  Vt/L,  after  which  the  initial  square  wave  (2.1)  of  width  20  . 
cells  has  been  carried  across  160  cells.  The  curves  (al) ,  (a2) ,  (a3) , 
correspond  to  the  schemes  of  Godunov  (which  coincides  in  the  linear 

*r. 

case  with  one-sided  differencing),  MacCormack  (equivalent  to  the  scheme 

«  ' 

of  Lax  and  Wendroff),  and  Eremin-Lipnitskii;  (b)  to  the  SHASTX  scheme;  (c). 


(d),  (e)  correspond  to  the  monotonie  MacCormack  scheme  with  smoothing 
according  to  algorithms  (1.7),  (1.10),  and  (1.9).  The  best  results 
were  obtained  with  SHASTX.  The  SHASTA^  algorithm  and  the  monotonic 
MacCormack  scheme  with  the  smoothing  (1.9)  were  very  slightly  inferior 
to  this.  The  smoothing  procedures  (1.7)  and  (1.10)  give  worse  results, 
which  are  comparable,  however,  with  the  results  obtained  from  the  third- 
order  scheme.  The  first-order  and  second-order  schemes  yielded 
unsatisfactory  approximations  to  the  solution. 

Since  the  solutions  of  the  linear  advection  equation  using  the 
algorithms  (1.7)  and  (1.10)  were  inferior,  however  slightly,  to  the  algorithm 
(1.9),  in  the  remainder  of  this  section  we  will  present  the  results  of 
calculations  carried  out  using  the  latter.  He  must  point  out,  however, 
that  these  other  algorithms  will  be  analyzed  below,  inasmuch  as  they 
possess  definite  advantages  in  solving  real  gas-dynamic  problems.  The 
results  of  the  present  section  are  applicable  in  problems  where  linear 
equations  are  solved,  for  example,  in  certain  meteorological  problems. 

Figure  2  illustrates  the  "flexibility”  of  the  various  schemes  in  the 
sense  of  Ref.  2,  i.e.,  the  capability  of  solving  problems  in  which  the 
Courant  number  varies  greatly  over  the  computational  region.  Here 
profiles  are  shown  at  t  *  3,  obtained  for  Courant  number  v*  0.6 
(curve  1)  and  v=  0.2  (curve  2)  using  first-order  (a),  second-order  (b) , 
third -order  (c)  and  the  monotonic  MacCormack  (d)  schemes.  The 
superiority  of  the  latter  is  indubitable.  The  calculation  was  not 
carried  out  with  SHASTX,  which  requires  that  condition  v<  0.5  be 
satisfied^.  « 

As  we  remarked  in  Section  1,  the  smoothing  which  removes  numerical 


oscillations  smears  out  physical  extrema.  Consequently,  the  test 
problem  using  initial  conditions  in  the  form  of  a  triangular  profile 
(2.2)  constitutes  a  stringent  test  of  a  scheme.  Figure  3  shows 
profiles  at  time  t  ■  1  (curves  1)  and  t  ■  20  (curves  2)  for  the  same 
schemes  as  in  Fig.  2,  if  v-  0.2.  The  first-order  scheme  practically 
"forgets"  the  Initial  conditions;  the  second-order  scheme  qualitatively 
changes  the  shape  of  the  pulse  as  time  elapses,  while  the  thid -order 
scheme  preserves  the  "height"  of  the  dlstrubance  well.  The  monotonic 
MacCormack  scheme  is  only  negligibly  inferior  to  SHASTX  (the  dashed  lines 


3.  The  evolution  of  a  gas-dynamic  shock 

The  problem  of  the  evolution  of  a  gas  dynamic  shock  involves  all  the 
characteristic  features  of  one-dimenrional  ideal  gas  flow:  creation  and 
propagation  of  a  shock  wave,  contact  discontinuity  and  rarefaction  fan. 

By  comparing  the  results  of  the  calculation  with  the  analytic  solution, 
ve  can  check  both  the  accuracy  with  which  the  characteristic  flow  regions 
are  transported  and  the  speed  with  which  they  develop-.  A  comparison 

of  the  various  schemes  mentioned  above  will  be  presented,  based  on  this 

\ 


problem. 

As  was  stated  earlier,  unphysicaloscillations  in  the  neighborhood  of 
shocks  and  steep  gradients  appear  because  of  the  nonvanishing  dispersion 


associated  with  a  difference  scheme,  the  nonlinearity  of  the  gas-dynamic 

24 


equations,  and  also  for  several  other  reasons  .  Besides  this,  it  is 
possible  that  the  various  nonlinear  smoothing  operators  may  Interact 

with  one  another  nonself consistently. 

The  system  of  equations  has  the  following  form: 


r*2  r*  2  •*' 

f  p|  d*+  f  pul' 
't,  Jt,  'x. 


dt  =  0, 


f2  2  f2 

\  C1  \  Xl 


dt  -  0. 


rx2  ,c2  rl2  p 

J  e  I  dx  +  I  (e  +  p  )u  I 


dt  *  0. 


Here  e  *  p/(y-l)  +  pu  /2,  where  y  is  the  adiabatic  index.  We  consider 

<- 

the  evolution  of  a  shock  with  initial  pressure  ratio  p^/p^  “  2  and 
density  ratio  pj/p^  "  and  Y  “  1.4.  In  Fig.  4  are  shown  density 
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profiles  obtained  using  the  Godunov  (a),  Eremin-Lipnitskii  (b), 

MacCormack  (c)  and  SHASTX  (d)  schemes.  The  exact  solution  is 
indicated  by  the  heavy  continuous  line  and  the  numerical  results 
up  to  t  B  25  by  a  fine  continuous  line.  The  first-order  scheme,  as 
one  might  expect,  badly  smears  the  shock  front  and  contact  discontinuity 
and  represents  the  rarefaction  fan  poorly.  The  Eremin-Lipnitskii  and 
MacCormack  schemes  are  nonmonotonic  and  yield  oscillations  in  the 
neighborhood  of  the  shocks.  SHASTX,  which  represents  the  shock  front 
very  well,  gives  rise  to  small  undershoots  in  the  density  in  the  neighborhood 
of  the  contact  discontinuity  and  rarefaction  fan.  For  this  case  the 
SHASTA  scheme,  as  described  by  Boris  and  Book,  gives  practically  the 
same  results  as  does  SHASTX.  All  of  the  schemes  accurately  represent 
the  pressure  profile  between  the  shock  wave  and  the  rarefaction  fan. 

In  Fig.  5  are  shown  results  of  calculations  using  the  monotonic 
MacCormack  scheme  with  various  forms  of  smoothing.  Fig.  5a  corresponds 
to  algorithm  (1.7).  One  shortcoming  of  this  approach  lies  in  the  presence 
of  undershoots  in  density  in  the  neighborhood  of  the  shocks  and  rarefaction 
fans.  Very  similar  results  are  obtained  when  one  applies  the  smoothing 
(1.9)  to  all  of  the  equations  in  system  (3.1)  (Fig.  5b). 

In  the  present  work  it  is  proposed  to  carry  out  self-consistent 
smoothing  of  the  density,  momentum,  and  energy  in  regions  where 
nonmonotonicity  in  the  density  profile  exists,  using  the  procedures 
(1.9)  and  (1.10)  described  above.  The  justification  for  this  lies  in 
the  fact  that  every  stable  one- dimensional  shock  has  a  density 
discontinuity.  In  particular,  this  self-consistent  smoothing 
introduces  no  changes  in  the  velocity  at  a  contact  discontinuity,  since 
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the  density  and  radnentum  are  smoothed  in  a  "related"  manner.  It  also 
somwehat  reduces  the  expenditure  of  machine  time,  since  in  order  to 
determine  nonmonotonicity,  one  array  of  numbers  instead  of  three  is 
searched  [Cf.(l.ll) j.  We  note  that  this  type  of  smoothing  cannot  be 
implimented  using  algorithms  (1.4)  -  (1.8)  because  of  their  two-stage 
character . 

The  density  profiles  obtained  using  self-consistent  smoothing  at 

til  / 

the  n  and  (n+1)  st  levels  are  shown  in  Figs.  5c  and  5d.  These 

results  are  clearly  superior  to  all  those  shown  previously.  It  is 

interesting  to  observe  that  the  difference  in  the  results  obtained  by 

T 

smoothing  according  to  the  algorithms  (1.9)  and  1.10)  is  much  smaller 
than  in  the  case  of  linear  advection  (Fig.  1).  This  is  indicative  of 
the  decisive  role  played  by  nonlinear  effects  in  numerical  solutions 
of  the  gas-dynamic  equations. 

We  pause  to  draw  attention  to  some  peculiarities  in  the  solution? 
we  have  found.  Most  of  the  calculations  were  run  with  At  «*  0.5. 

SHASTA,  SHASTX,  and  the  MacCormack  scheme  with  self-consistent 
smoothing,  however,  yielded  a  highly  oscillatory  solution.  The 
oscillations  are  a  consequence  of  the  nonself consistency  of  the  smoothing 
and  the  interaction  between  the  resulting  errors  in  the  various  fluid 
quantities.  A  separate  paper  will  be.  devoted  to  a  detailed  analysis 
of  the  nonlinear  properties  of  the  different  schemes,  the  nonmonotonic! ty 
resulting  from  various  causes,  and  ways  of  contending  with  the  unphysical 
oscillations  that  result.  We  therefore  will  not  go  inot  the  reasons  for 
the  appearance  of  computational  ripples  in  detail,  but  content  ourselves 
with  remarking  that  they  vanish  as  the  time  step  is  reduced.  The  > 
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solutions  shewn  in  Fig-  4c  and  5b  were  obtained  with  a  time  step 
t  -  0.25. 

We  now  compare  tin;  respective  running  times.  Taking  as  the  unit  of 
time  that  required  for  a  single  step  according  to  the  HacCormack  scheme, 
the  corresponding  running  times  for  the  other  schemes  ate  as  follows: 
monotonic  HacCormack  with  self  consistent  smoothing,  1.4;  with 
nonselfconsistent  smoothing  according  to  algorithm  (1.9),  1.5;  using 
algorithm  (1.7),  2;  for  the  Godunov  scheme  2  -  2.5,  depending  on  how 
the  "large"  quantities  are  calculated;  for  SHASTA  and  SHAoTX,  3;  and 
for  the  Eremln-Lipnitskii  scheme,  2.5. 

We  next  consider  an  example  in  which  a  strong  shock  wave  is  formed. 
The  initial  parameters  are  taken  to  he  the  same  as  in  the  paper  of  Boris 
and  Bookli(p2/p1  »  480,  p^/p^  **  8,  y  «  5/3,  At:  «  0.02).  The  Eremin- 
Lipnitskii  and  HacCormack  schemes  in  the  absence,  of  an  additional 
artificial  viscosity  do  not  permit  calculation  of  such  a  flow.  SDA'-T/. 
and  the  monotonic  MacCormack  scheme  with  smoothing  according  to 
algorithm  (1.7)  give  approximately  the  seme  results  as  SHaSTX .  We 
therefore  make  further  comparisons  of  the  results  obtained,  using  the 
method  of  Godunov,  SHASTX,  and  me no tonic  MacCormack  scheme  with  h>  ’■ 
consistent  smoothing  given  by  (1.9),  In  Figs.  6  and  7  ai  c  shown  eensih* 
and  pressure  profiles  obtained  using  these  r.ohemor,  at  t.im-  t  **  2  (deshrd 
lines  and  t  =  4  (fine  continuous  lines).  A  vertical  dashed*  dotted  lint 
is  used  to  show  the  initial  position  of  the  shock.  The  virtues  and 
shortcomings  of  the  numerical  methods  under  invest igat  ton  arp 
clearly  shown  i..  Table  1,  in  which  the  values  of  density  and  relative 
errors  at  time  t  »  4  are  shown.  In  the  left  column  of  the  table  is  shown 


the  number  of  computational  cells.  The  Godunov  scheme  strongly 
smears  the  shock  front  and  contact  discontinuity  and  approximates  the 
exact  solution  poorly  In  the  region  of  the  rarefaction  fan.  SHASTX 
is  somewhat  better  than  the  monotonic  MacCormack  scheme  in  representing 
the  shock  wave  and  rarefaction  fan,  but  somewhat  inferior  to  it  in 
representing  the  contact  discontinuity.  In  addition,  SHASTX  gives 
incorrect  values  in  the  neithborhood  of  the  initial  shock. 

The  monotonic  HacCormack  scheme  and  SHASTX  have  good  dynamical 
properties.  They  permit  relatively  fast  determination  of  the  true  flow 
pattern  (Cf.  the  Godunov  scheme).  This  last  property  is  especially  important 
in  solving  nonstationary  problems. 


SI 
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V 


A.  Calculation  of  three  dimensional  interaction  of  a  shock  wave  with 


an  obstacle. 

The  possibilities  of  the  proposed  method  are  illustrated  below  for 

the  example  of  reflection  of  a  plane  shock  wave  from  a  obstacle  of 

finite  width.  The  initial  stage  of  the  interaction  precess  is 

investigated,  in  which  the  reflected  shock  wave  begins  to  form  and  the 

pressure  loading  on  the  object  is  a  maximum.  We  remark  that  this  problem 

is  of  interest  in  its  o\m  right,  inasmuch  as  experimental  determination  of 

nonstationary  three-dimensional  flows  encounters  serious  difficulties. 

The  problem  is  solved  in  Cartesian  coordinates.  The  mesh  is 

uniform.  Timestep  splitting  was  employed  in  solving  the  difference 
25  26 

equations  *  which  permitted  a  substantial  saving  in  running  time. 

The  conservation  laws  in  Eulerian  coordinates  have  the  form 


^"pj2dv+  J  2  (j)  pV*n 

V  C1  tl  s 

/s5l  2->-  */  i  (pW  +  p) 

V  *1  fcl  s 


r  i"2  A  r  %  -  - 

J  e  I  dv  +  I  (p  (e  +  p)  V  •  n 

1  t,  s 


ds  dt  *=  0, 


n  ds  dt  =  0, 


ds  dt  »•  0 


We  denote  by  (At)  the  operator  which  performs  a  step  in  the  direction  o: 


where 


ruk  -  ‘i"k  -  £  ^s-Fuk), 

n+l/3  b  1  r.  n  ~  ft  _  f  \-j 

‘ijk  2  L  ijk  ijk  Ah  K‘ ijk  i-q,j-r,k-s' J» 


pu  u  +  t  p, 
y  a 

pu'  s 


q“5,r“6,s«»6 

xa  ya  za. 

The  smoothing  is  carried  out  according  to  the  algorithm  (1.9).  The 
transverse  momentum  components  were  smoothed  independently. 

The  complete  transport  operator  will  look  as  follows: 

L(2At)  «  L  (At)  L  (At)  L  (At)  L  (At)  L  (At), 
y  x  z  x  / 

Written  thus,  even  if  the  split  operators  LQ(At)  do  not  commute,  the 
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difference  scheme  has  second-order  accuracy  over  a  double  timestep 
The  stability  condition  for  the  split  equations  is 


where 


At  £  min (At  ) ,  a  *  x,y,z, 
a  a 


At  S  min  (~ — 

°  ijk  ^ 


and  a  is  the  speed  of  sound.  This  condition  turns  out  fco  be 
considerably  less  stringent  than  for  the  unsplit  operators: 


At  £  min  ( 


5 — 7=*- 

+  a/  3 


_1 


Since  At^  (  x  is  the  direction  of  propagation  of  the  shock  wave) 

is  usually  much  smaller  than  At  and  At  ,  it  makes  sense  to  carry  out 

y  z 

the  calculation  in  the  y  and  z  directions  with  a  double  timestep.  Thus 
using 

LQ(2At)  *  La<At)  Lft(At),  a  =  y,z. 


we  finally  obtain 


n+2mAt 


ijk 


L  (At)[L  (At)  L  (2At)  L  (At)  L  (2At)J 

y  x  z  x  y 


M-l 


.  L  (At)L  (2At)L  (At)L  (At)f  . 

x  z  x  y  X  J  K 


Use  of  this  algorithm  instead  of 


n+2mAt 


-  Lm(2At)f . . 


ijk 


ijk 


permits  reduction  of  the  computational  volume  by  a  factor  of  1.5. 

In  the  test  the  computationalmesh  consisted  of  20  X  20  X  15  cells. 

Calculation  of  a  single  case  up  to  t  =  1.2  -  1.5  required  about  2  hours  of 

time  on  the  M4030.  Here  t  =  tV  /h,  where  V  .is  the  speed  of  the 

amp  amp  r 

reflected  shock  wave  which  arose  as  a  result  of  the  interaction  of  the 
incident  wave  with  the  square  object,,  and  h  is  the  step  height. 

A  typical  flow  pattern  is  shown  in  Fig.  8.  At  t  =  1.1  the  isobars 
are  shown  in  (a)  with  separation  Ap  =  10.0, and  the  lines  of  constant 
Mach  number  (b)  at  intervals  AM  =  0.2.  The  Mach  number  of  the  incident 
shock  wave  M  =5.0  and  the  ratio  of  width  to  height  of  the  step  was 
1/ h  =  5.2.  The  adiabatic  index  was  y  =  1.4  and  tha  pressure  in  the 
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V 


undisturbed  gas  was  =  1.0.  It  is  clear  that  the  method  permits 
excellent  resolution  of  the  shock  wave  and  the  flow  in  the  neighborhood 
of  the  corner  points. 

In  Fig.  9  are  shown  the  isobars  in  two  cross  sections  x  *  const, 
located  equidistantly  and  parallel  to  the  boundary  of  the  step  at  a 
distance  0.3  h.  The  ratio  I/h  =■  2.0  in  (a),  3.6  in  (b) ,  and  5.2  in  (c). 

In  Fig.  9a  are  shown  the  isobars  for  the  plane  case  (i/h  «  ») . 

Figure  10  illustrates  the  influence  of  the  ratio  i/h  on  the  bending 
of  the  shock  wave  in  the  xy  plane.  Curves  1,2,3  ^correspond  respectively 
to  l/h  =  2/0,  3.6,  and  5.2  (the  bending  is  related  to  the  bending  for  the 
plane  case) . 

The  authors  thank  Yu.  P.  Golovachev  for  his  unstinting  attention  to 
this  work  and  useful  discussions,  and  also  V.M.  Goloviznin,  G.L.  Stenc.hlknv, 
and  A.P.  Favorskii  for  valuable  suggestions. 
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a  (Godunov  nethod) 


b  (SHASTX) 


c  (nonotonic 
MacCormack  scheme) 
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for  Intended  Recipient) 

Olcy  Attn  Mail  Ser  Sec  W  Roherty 

Olcy  Attn  Mail  Ser  Sec  3141 

Olcy  Attn  Mail  Ser  Sec  L  Vortman 

Olcy  Attn  Mail  Ser  Sec  A  Chaban 

Olcy  Attn  Mail  Ser  Sec  L  Hill 

Science  Applications,  Inc. 

P  0  Box  2351 
La  Jolla,  CA  92038 

Olcy  Attn  Technical  Library 

Science  Applications,  Inc. 

1257  Tasman  Drive 
Sunnyvale,  CA  94086 

Olcy  Attn  J  Dishon 

Science  Applications,  Inc. 

2450  Washington  Avenue 
San  Leabdro,  CA  94577 
Olcy  Attn  D  Maxwell 

Olcy  Attn  D  Bernstein 

Science  Applications,  Inc. 

P  0  Box  1303 
McLean,  VA  22102 

Olcy  Attn  M  Knasel 

Olcy  Attn  B  Chambers  III 

Olcy  Attn  R  Sievers 

Olcy  Attn  J  Cockayne 

Southwest  Research  Institute 
P  0  Drawer  28510 
San  Antonio,  TX  78284 
Olcy  Attn  W  Baker 

Olcy  Attn  A  Wenzel 

SRI  International 
333  Ravenswood  Avenue 
Menlo  Park,  CA  94025 

Olcy  Attn  G  Abrahamson 

Systems,  Science  £.  Software,  Inc. 

P  0  Box  1620 
La  Jolla,  CA  92038 

Olcy  Attn  Library 
Olcy  Attn  D  Grine 
Olcy  Attn  T  Riney 
Olcy  Attn  K  Pyatt 


Teledyne  Brown  Engineering 
Cummings  Research  Park 
Huntsville,  AL  35807 

Olcy  Attn  J  Ravenscraft 

Terra  Tek,  Inc. 

420  Wakara  Way 
Salt  Lake  City,  UT  84108 
Olcy  Attn  Library 
Olcy  Attn  S  Green 
Olcy  Attn  A  Jones 


Tetra  Tech,  Inc. 

630  H  P.osemead  Blvd. 
Pasadena,  CA  91107 
Olcy  Attn  L  Hwang 
Olcy  Attn  Library 


TRW  Defense  &  Space  Sys  Group 
One  Space  Park 
Redondo  Beach,  CA  90278 


01c 

Clc 

02cv 

Olcy 

01c 

01c 


Attn  I  Alber 
Attn  Tech  Infor  Ctr 
Attn  N  Lipner 
Attn  P  Bhutta 
Attn  D  Baer 
Attn  R  Plebuch 


TRW  Defense  6  Space  Sys  Group 
P  0  Bo:  1310 

San  Bernardino,  CA  92402 
Olcy  Attn  E  Wong 
Olcy  Attn  P  Dai 


Universal  Analytics,  Inc. 
7740  W  Manchester  Blvd 
Play a  Del  Rev,  CA  90291 
Olcv  Attn  E  Field 


Weidlinger  Assoc.,  Consulting  Eng 
110  E  59th  Street 
Sew  York,  NY  10022 

Olcy  Attn  M  Baron 

Weidlinger  Assoc.,  Consulting  Eng 
3000  Sand  Hill  Road 
Menlo  Park,  CA  94025 

Olcy  Attn  J  Isenberg 

Westinghouse  Electric  Corp. 

Miarine  Division 
Hendy  Avenue 
Sunnyvale,  CA  94088 
Olcy  Attn  W  Volz 


44 


